/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#include <math.h>
#include <stdlib.h>
#include <string.h>

#include "fastonebigheader.h"

#define exp fastexp
#define log fastlog

#ifdef __SSE2__
static float
denominator_acc_loop (float*       logp,
                      float        max,
                      unsigned int n)
{
  float sum = 0.0;
  v4sf vsum = _mm_setzero_ps ();
  v4sf vmax = _mm_load1_ps (&max);
  unsigned int i;

  for (i = 0; i < n && (((uintptr_t) logp) % 16) > 0; ++i)
    {
      sum += fastexp (logp[i] - max);
    }

  for (; i + 4 < n; i += 4)
    {
      v4sf varg = _mm_load_ps (logp + i);
      vsum += vfastexp (varg - vmax);
    }

  for (; i < n; ++i)
    {
      sum += fastexp (logp[i] - max);
    }

  return sum + v4sf_index (vsum, 0) +
               v4sf_index (vsum, 1) +
               v4sf_index (vsum, 2) +
               v4sf_index (vsum, 3);
}
#else
#warning "not using SSE2"
static float
denominator_acc_loop (float*       logp,
                      float        max,
                      unsigned int n)
{
  float sum = 0.0;

  unsigned int i;

  for (i = 0; i < n; ++i)
    {
      sum += fastexp (logp[i] - max);
    }

  return sum;
}
#endif

/* Ugh ... I should have used C++ */

#define logprating(_rating, _loga, _rank, _logbeta, _c, _d, _n_labels)  \
  ({                                                                    \
    /*                                                                  
     * ok we want to compute the logarithm of
     * 
     * exp (exponents[rating])     exp (exponents[rating] - max)
     * ----------------------- = --------------------------------
     * sum_i exp (exponents[i])   sum_i exp (exponents[i] - max)
     *
     * the RHS is more numerically stable
     */                                                                 \
                                                                        \
    float _a[_rank];                                                    \
    for (unsigned int _logprn = 0; _logprn < _rank; ++_logprn)          \
      {                                                                 \
        _a[_logprn] = exp (_loga (_logprn));                            \
      }                                                                 \
    float _beta = exp (_logbeta);                                       \
    float _logp[_n_labels];                                             \
    unsigned int _maxindex = 0;                                         \
                                                                        \
    for (unsigned int _logpri = 0; _logpri < _n_labels; ++_logpri)      \
      {                                                                 \
        _logp[_logpri] = 0;                                             \
        for (unsigned int _logprn = 0; _logprn < _rank; ++_logprn)      \
          {                                                             \
            _logp[_logpri] -= _a[_logprn] *                             \
                              _beta *                                   \
                              _c (_logprn) *                            \
                              _d (_logprn, _logpri);                    \
          }                                                             \
                                                                        \
        if (_logp[_logpri] > _logp[_maxindex])                          \
          {                                                             \
            _maxindex = _logpri;                                        \
          }                                                             \
      }                                                                 \
                                                                        \
    float _max = _logp[_maxindex];                                      \
    float _denominator = 0;                                             \
                                                                        \
    _denominator = denominator_acc_loop (_logp, _max, _n_labels);       \
                                                                        \
    _logp[_rating] - _max - log (_denominator);                         \
  })

#define dlogpratingdall(_dloga, _dc, _dd, _rating,                      \
                        _loga, _rank, _logbeta, _c, _d, _n_labels)      \
  ({                                                                    \
    float _a[_rank];                                                    \
    for (unsigned int _logprn = 0; _logprn < _rank; ++_logprn)          \
      {                                                                 \
        _a[_logprn] = exp (_loga (_logprn));                            \
      }                                                                 \
                                                                        \
    float _beta = exp (_logbeta);                                       \
    float _logp[_n_labels];                                             \
    unsigned int _maxindex = 0;                                         \
                                                                        \
    for (unsigned int _logpri = 0; _logpri < _n_labels; ++_logpri)      \
      {                                                                 \
        _logp[_logpri] = 0;                                             \
        for (unsigned int _logprn = 0; _logprn < _rank; ++_logprn)      \
          {                                                             \
            _logp[_logpri] -= _a[_logprn] *                             \
                              _beta *                                   \
                              _c (_logprn) *                            \
                              _d (_logprn, _logpri);                    \
          }                                                             \
                                                                        \
        if (_logp[_logpri] > _logp[_maxindex])                          \
          {                                                             \
            _maxindex = _logpri;                                        \
          }                                                             \
      }                                                                 \
                                                                        \
    float _max = _logp[_maxindex];                                      \
    float _denominator = 0;                                             \
    float _ddenominator_dloga[_rank][_n_labels];                        \
    float _ddenominator_dlogbeta = 0;                                   \
    float _ddenominator_dc[_rank][_n_labels];                           \
    float _ddenominator_dd[_rank][_n_labels];                           \
                                                                        \
    memset (_ddenominator_dloga,                                        \
            0,                                                          \
            sizeof (_ddenominator_dloga));                              \
                                                                        \
    memset (_ddenominator_dc,                                           \
            0,                                                          \
            sizeof (_ddenominator_dc));                                 \
                                                                        \
    memset (_ddenominator_dd,                                           \
            0,                                                          \
            sizeof (_ddenominator_dd));                                 \
                                                                        \
    for (unsigned int _logpri = 0; _logpri < _n_labels; ++_logpri)      \
      {                                                                 \
        float _d = _logp[_logpri] - _max;                               \
        float _expd = exp (_d);                                         \
        _denominator += _expd;                                          \
        _ddenominator_dlogbeta += _expd * _d;                           \
                                                                        \
        for (unsigned int _logprn = 0; _logprn < _rank; ++_logprn)      \
          {                                                             \
            _ddenominator_dloga[_logprn][_logpri] -=                    \
                                            _expd *                     \
                                            _a[_logprn] *               \
                                            _beta *                     \
                                            _c (_logprn) *              \
                                            _d (_logprn, _logpri);      \
            _ddenominator_dloga[_logprn][_maxindex] +=                  \
                                            _expd *                     \
                                            _a[_logprn] *               \
                                            _beta *                     \
                                            _c (_logprn) *              \
                                            _d (_logprn, _maxindex);    \
            _ddenominator_dc[_logprn][_logpri] -=                       \
                                            _expd *                     \
                                            _a[_logprn] *               \
                                            _beta *                     \
                                            _d (_logprn, _logpri);      \
            _ddenominator_dc[_logprn][_maxindex] +=                     \
                                            _expd *                     \
                                            _a[_logprn] *               \
                                            _beta *                     \
                                            _d (_logprn, _maxindex);    \
            _ddenominator_dd[_logprn][_logpri] -=                       \
                                            _expd *                     \
                                            _a[_logprn] *               \
                                            _beta *                     \
                                            _c (_logprn);               \
            _ddenominator_dd[_logprn][_maxindex] +=                     \
                                            _expd *                     \
                                            _a[_logprn] *               \
                                            _beta *                     \
                                            _c (_logprn);               \
          }                                                             \
      }                                                                 \
                                                                        \
    memset (_dloga, 0, _rank * sizeof (_dloga[0]));                     \
    memset (_dc, 0, _rank * sizeof (_dc[0]));                           \
    for (unsigned int _logpri = 0; _logpri < _n_labels; ++_logpri)      \
      {                                                                 \
        for (unsigned int _logprn = 0; _logprn < _rank; ++_logprn)      \
          {                                                             \
            _dloga[_logprn] -=                                          \
              _ddenominator_dloga[_logprn][_logpri] / _denominator;     \
            _dc[_logprn] -=                                             \
              _ddenominator_dc[_logprn][_logpri] / _denominator;        \
            _dd[_logprn * _n_labels + _logpri] =                        \
              -_ddenominator_dd[_logprn][_logpri] / _denominator;       \
          }                                                             \
      }                                                                 \
                                                                        \
    for (unsigned int _logprn = 0; _logprn < _rank; ++_logprn)          \
      {                                                                 \
        _dloga[_logprn] -= _a[_logprn] *                                \
                           _beta *                                      \
                           _c (_logprn) *                               \
                           _d (_logprn, _rating);                       \
        _dloga[_logprn] += _a[_logprn] *                                \
                           _beta *                                      \
                           _c (_logprn) *                               \
                           _d (_logprn, _maxindex);                     \
        _dc[_logprn] -= _a[_logprn] *                                   \
                        _beta *                                         \
                        _d (_logprn, _rating);                          \
        _dc[_logprn] += _a[_logprn] *                                   \
                        _beta *                                         \
                        _d (_logprn, _maxindex);                        \
        _dd[_logprn * _n_labels + _rating] -=                           \
                                 _a[_logprn] *                          \
                                 _beta *                                \
                                 _c (_logprn);                          \
        _dd[_logprn * _n_labels + _maxindex] +=                         \
                                   _a[_logprn] *                        \
                                   _beta *                              \
                                   _c (_logprn);                        \
      }                                                                 \
                                                                        \
    _logp[_rating] - _max - _ddenominator_dlogbeta / _denominator;      \
  })

